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Abstract: We present a novel method to implement time-delayed propagation of radiation fields in cosmo¬ 
logical radiative transfer simulations. Time-delayed propagation of radiation fields requires construction 
of retarded-time fields by tracking the location and lifetime of radiation sources along the corresponding 
light-cones. Cosmological radiative transfer simulations have, until now, ignored this “light-cone effect” 
or implemented ray-tracing methods that are computationally demanding. We show that radiative trans¬ 
fer calculation of the time-delayed fields can be easily achieved in numerical simulations when periodic 
boundary conditions are used, by calculating the time-discretized retarded-time Green’s function using the 
Fast Fourier Transform (FFT) method and convolving it with the source distribution. We also present a 
direct application of this method to the long-range radiation field of Lyman-Werner band photons, which 
is important in the high-redshift astrophysics with first stars. 
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1. Introduction 

Cosmological radiative transfer simulations are used to 
study astrophysical processes which occur on a large, 
cosmic scale. A notable example of such processes 
is the process of cosmic reionization, in which indi¬ 
vidual H II regions take up large volumes, of order 
(~ 20 comoving Mpc)^ or larger, after ~ 50% of all 
the bar yons are reionized by astrophysical radiation 
sources (iFurlanetto fc Ohll2005l: iFiirlanetto et al.l[200^ 

llliev et al.ll2007t IZahn et al. l2007ll ! The typical radi¬ 

ation sources long after the recombination epoch are 
stars and quasars, which emit predominantly ultra¬ 
violet (UV) and X-ray photons, respectively. The mean 
free path of hydrogen-ionizing (H-ionizing) UV pho¬ 
tons is simply the average size of H II regions^ while 
the mean free path of X-ray photons is much larger 
than that of UV photons due to the relatively small 
optical depth (the scattering cross section of UV pho- 
tons is much larger t han that of X- ray photons: e.g. 
iMesinger et al.l [20131 : IXu et al.ll2014ll . Therefore, the 
mean free path of X-ray photons is truly cosmological, 
surpassing several hundred Mpc whe n the rest-fram e 
photon energy is larger than a few keV (|Xu et al.ll2014ll . 
Another example of long mean free paths is the Lyman- 
Werner (LW) band photons. These photons lie at the 
energy range of 11 - 13.6 eV, and thus traverse cos¬ 
mological distances freely until they redshift into hy¬ 
drogen Lyman resonance lines and are scattered off 
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■'■Note that even after the intergalactic medium (IGM) is fully 
ionized, which is believed to occur at z < 7, the mean free path 
of UV photons is severely restricted to a few comoving Mpc. 
This is due to the rapid increase of the number of Lyman limit 
systems, which contain neutral hydrogen atoms inside and thus 
consume hydrogen-ionizing UV photons (e.g. IMcOuinn et all 
120111 : 1 Alvarez fc Abeill20i;^ . 


neutr al hydrogen atoms (jHaiman et aLl2d00tlAhn et al.l 
I 2 OO 9 II . Practically all LW band photons are scattered 
multiple times and reprocessed into lower-energy pho¬ 
tons when they traverse the “Lyman-Werner horizon” 
flw ~ 100 [21/(1 -I- Mpc, w here z is the reds hift 
at which the photons are emitted ( Ahn et ^l2009t) . 

The light-cone effect, which denotes the effect related 
to fields travelling at the speed of light and thus tracing 
past-time events along the corresponding light cones, 
is inherent in any physical processes on cosmic scales. 
Radiation fields and gravitational fields are the obvious 
examples governed by the lig ht-cone effect: both fields 
travel at the speed of light (see lHwang et al.l[2008l show- 
ing that the gravitational field travels at the speed of 
light in the form of the electric part of the Weyl ten¬ 
sor). The relevant time or length scale over which the 
light-cone effect matters can be estimated by an effec¬ 
tive time-scale tic = //(df/dt), where / is any physical 
quantity relevant to the problem of interest and df/dt 
is the change rate of /. For astrophysical problems, 
the time scale may be naively estimated by the life¬ 
time of radiation sources of interest. In contrast, grav¬ 
itational fields under non-relativistic motion (velocity 
<C c) will be corrected from the Newtonian field, which 
is constructed in an action-at-a-distance manner, only 
very slightly (about 10“® - 10“^ or somewhat larger if 
cumulative effects are considered) by gen eral relativis- 
tic effects including this light-cone effect (iHwang et al.l 
1200811 . 


Light-cone effect of radiation fields traversing cos¬ 
mological distances has been usually approximated 
by a uniform, spatially- averaged quantity in semi- 
analytical studies (e.g. iHaiman et al.l l200fll for LW 
and H-ionizing backgrounds), or calculated by com¬ 
putationally expensive, brute-force met hods of ray- 
tracing with finite speed of light (e.g. iBrvan et al.l 
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I 2 OI 4 I for H-ionizing background). The first account 
of the inhomogeneous LW backgroun d in n umerical 
simulations was taken by lAhn et ^ ( 20091! . whose 
consistency with semi-analytical and semi-numerical 
approach was confirmed elsewhe re (iDiikstra et al.l 


20081: Holzbauer fc Furlanettol l2012l! . The method of 
Ahn et al.l (|2009l! was later used in simulating the pro¬ 


cess of cos mic reionization under the influence of the 
first stars ( Ahn et ^120121: iFialkov et al.ll201,'l^ . This 
method has been improved for faster computation using 
FFT. Its application to other radiation fields affected 
by the light-cone effect is straightforward, as long as 
the field strength from the source is only a function of 
distance but not of direction. Even though this im- 
proved me t hod w as already used in the simulation by 
lAhn et ^ ( 20121! . its formalism and generality have 
never been described explicitly. Thus we lay out the 
generic formalism, a numerical scheme and the actual 
application i n this paper. A similar account for X-ray 
ba ckground (Xu et al.ll20ll! used the same machinery 
of lAhn et al. (l2009l! with modifications to implement 
the opacity due to hydrogen and helium atoms. 

This paper is organized as follows. In Section [5J 
we describe the generic formalism and the numerical 
scheme. In Section [31 we show how the LW background 
is calculated using this scheme, and also quantify the 
possible error from ignoring the light-cone effect. In 
Section m we conclude our findings and discuss further 
prospects and limitations. 


2. Retarded Radiation Field: Formalism 
AND Numerical Method 

2.1. Green’s Function Formalism 

A generic radiation field E(x, t) at position x and time 
t in a homogeneous and isotropic background spacetime 
satisfies a d’Alembertian-type equation. In the trivial 
case of a Minkowski back ground, electromagnetic (EM) 
waves (see l,Iacksonlll998l! follow 

t) =-47r/(x, t), (1) 

where F'(x, t) is either the scalar potential $ or the 
vector potential A, and /(x, t) is the corresponding 
source term (the charge density or the current density, 
respectively). In this case, F'(x, t) can be cast into an 
convoluted integral form 

F(x, t) = y (fx'dt'G{x, t; x', t'), (2) 

where the retarded Green’s functiorQ (RGF henceforth) 
C?(x, t; x', t') satisfies 

t; x', t') = -47rd(x - x.')S{t - t') 

"" _W 

^We take the simple stance that the initial condition is well es¬ 
tablished in astrophysical problems, and therefore the retarded 
Green’s function is more adequate than the advanced Green’s 
function. 


and is given in a closed form as 


G(x, t\ x', t') 



( 4 ) 


Note that Equation m is in the most generic form. 
The light-cone effect or causality is apparent in equa¬ 
tions @ and o, as the field at observing time t is 
caused by a source along the light cone at past time 
t' = t — |x —x'l/c. Even though equations ([T]), ([3]) 
and are applicable both to $ and A only in the 
Lorentz gauge, it is well known that the electromag¬ 
netic field, as a physical quantity derived from both $ 
and A, still conserve s this causality under any gauge 
choice ( ,Iacksonll2002l: see also a nice example of Prob¬ 
lem 6.20 in lJacksonlllQ^ . 


Let us consider a generic field E(x, r) in the space- 
time described by the Robertson-Walker metric with 
null curvature. 


ds^ = <^dt^ — a^{t)dx^ = a{T) {dr^ — dx^) , (5) 


where x is the comoving spatial coordinate, a{t) is the 
scale factor at conformal cosmic time r with normal¬ 
ization a = 1/(1 -I- z), and r is defined by 



It is convenient to use r and x because the null geodesic 
is simply given by d |x| /dT=l. One can then expect 
that the generic RGE will take the following form: 


G(x, r; x', r') = d (r'- [r - |x - x'|]) 

xJ'dx-x'l , T, r'), (7) 


where J'dx — x'|, r, r') is a generic, geometrical dilu¬ 
tion factor. Note that J^(|x —x'|, r, t'), in general, 
reflects cosmological effects such as photon-redshifting, 
and therefore depends also on r and r'. The corre¬ 
sponding integral form will become (c.f. Equation [3]) 


F(x, t) = J d^x' dr' r; x', r')/(x', r'), 


( 8 ) 


where /(x', r') is the source term properly defined in 
the (x, t) coordinate system. 

We now decompose Equation ([8]) into discrete terms 
which are affected by different time slices, and also 
apply a periodic boundary condition to assimilate the 
usual simulation set-up. Equation (|8]) can be rewritten 
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as 


F(x, r) 

= X! / i '^')/(x', t') 

= / d^cc' f dr '(5 (r'— [t — |x — x'l]) 

n 

xJ^(|x-x'|, T, t')/(x', t') 

T / d’T IF f 

„ d V ''■n-Fl - T-n J 

xJ^(|x-x'|, T, t')/(x', T„) 

= ^ / d^x'J'ndx-x'l; t)/(x', r„), 

n 

= ^F„(x, r) ( 9 ) 


where /(x', t„) approximates the source term during 
t' = [t„, t„+i] as a fixed 3D quantity, LE(a;) is a top- 
hat window function such that 


spatially convoluted integrals (Equation (TU]). In the 
next section, we will describe how we calculate the full 
retarded field under a uniform Eulerian grid. 

2.2. Numerical Method 

We now develop a numerical method to calculate the 
retarded fields. A scheme for calculating 3D fields in 
numerical simulations becomes the most transparent 
under a uniform grid, and thus we restrict the descrip¬ 
tion also to the case of a uniform grid. This choice is 
further justified by the fact that a retarded field origi¬ 
nates from sources located far from the observing point, 
whose impact on each observing point is greatly aver¬ 
aged out and thus can be accurately calculated even on 
a uniform grid. 

If we discretize relevant quantities on a uniform grid, 
Equation (US turns into 


Fn{i, j, k) 


N -1 N -1 N -1 


X! X! X! ^ ~ ^ ~ 

i'—O 0 k'—O 

Xfnil', f, k'), ( 11 ) 


W(x) 


1 if 0 < X < 1 
0 otherwise. 


and the RGF and the partial contribution to E(x, t) 
from the sources during r' = [r„, t„+i] are defined as 
J^„(|x —x'l; t) and F„(x, r), respectively. It is im¬ 
portant to choose small enough time gaps such that 
/(x', r') does not evolve much during each time step. 
Now, we apply the periodic boundary condition: the cu¬ 
bical domain with size L (volume V = L^) is infinitely 
repeated side to side and thus satisfies 


E„(x, r) 

nV 

= ^ / (fix'Tn[\x-{x' + {u, V, w)L}\-, t] 

U, V, W 

x/(x' -I- (u, V, w)L, T„) 

.V 

= / (fix' ^ Tn[\x-{x' + {u, V, w)L}\-, t] 

u, w 

x/(x', r„) 

pV 

= J d^x'Qnix-x')f{x', Tn), (10) 


where u, v, w are integer indices running from — oo to oo 
along the three axes of the cubical domain, both x and 
x' are confined to one domain, the second identity is 
guaranteed by the periodicity of the source term, and 
the last identity defines the “effective” RGF Gn{x. — 
x'). Note that for a given observing point (x, r) and 
the past time-slice at [t„, t„+i], the (u, w, w) points 
contributing to Gn(x. — x') are only those confined in 
the spherical shell bounded by comoving radii r — Tn+i 
and T — Tn- 

We have just cast the full integral (Equation [5]) for 
the retarded field into the summation (Equation iQ]) of 


where the domain is divided into cells with size 
I = L/N, {i, j, k} and {z', j', k'} are integer 
indices for real-space positions, and we absorbed 
fi into fn{i', j', k') by letting / d^x'f{x', t„) i— 
fnii', j', k'). If A{p, q, r) is the Fourier com¬ 
ponent of some field A(i, j, k) given by 


A(z, j, k) 


N-lN-lN-l 

E E E q, r) 

p—0 q—Q r=0 


X exp 


27ri(zp + jq + kr) 

N 


( 12 ) 


where the imaginary unit i is shown in bold face to 
avoid confusion with the index z, then the convolution 
theorem gives 


FnijP, q, r) = G„(p, q, r) fnijp, q, r), (13) 

where {p, q, r} are the indices of the Fourier-space co¬ 
ordinates. 

The merit of Equations (HH) - (USD is that the re¬ 
tarded field can now be calculated by EFT as long as 
periodicity is assumed. The algorithm then becomes 
straightforward, as follows: 

1. For given (observing) time r, divide the lookback 
time into discrete time slices with r' = [t„, t„+i]. 

2. For a given (past) time slice n, locate a unit source 
at the origin (z', j', k') = (0, 0, 0) of the real- 
space domain (box), and attach identical boxes 
- each containing one unit source at its origin - 
side by side. The total number of boxes to attach 
may be limited by the “effective horizon”, beyond 
which the contribution to Gnih J, k) (or to the net 
Fn{i, j, k) due to a rapid decrease of j', k') 
in lookback time) becomes too weak to include. 






4 


Ahn 



Figure 1. RGFs for LW intensity at observing redshift z = 27.90, in a periodic box of size 114//i (comoving) Mpc and a 
uniform grid of 256® cells. Both the real-space (top row) and Fourier-space (bottom row) RGFs {Qn and Qn respectively) are 
shown in log scales, on a slice containing 4 corners of the box. For Qn which has large dynamic range, we took the absolute 
value of its real part, |Re(5n)|, for better visualization. From left to right, RGFs of progressively longer lookback-time 
slices are shown, in units arbitrarily chosen but consistent throughout time slices. Lookback-time range is shown in terms of 
redshift, below each set of Qn and Qn- Each redshift range corresponds to 2 Myrs of cosmic time. Note that the right-most 
panels correspond to the time slice containing the horizon rnw. 


Calculate Gnii, j, k) from this configuration by a 
direct summation (see its definition in Equation 
[ini) over the contributions from these sources, when 
the functional form of J^(|x — x'|, r, r') in Equa¬ 
tion is known. Note that Gn{i, j, k) at each 
observing point {i, j, k} is affected only by the 
unit sources inside the spherical shell given by 
T — Tn +1 < ros < T — Tn, where Tos is the comoving 
distance between X = {i, j, k)l andx' = (u, v, w)L 
and the index notation follows those of Equation 
(fTUj) and (HD). 

3. Get Gn{p, <?, r) by applying EFT to Gnii, j, k) cal¬ 
culated from step 2. 

4. Get /(p, q, r) by applying EFT to /„(i, j, k). 

5. Get Fn{p, q, r) from Equation (TT^ . and then apply 
EFT to get j, k). 

6 . Finally, sum over time slices to get F{i, j, k) = 

j, k). 


be iVunit-source (fh/L)^, and the computation of each 
Gnii, j, k) will require ~ iVmesh(?')A^unit-source opera¬ 
tions. 


3. Application: Lyman-Werner Background 


As an application, we apply the method from Section 
to the calculation of LW background. LW bands 
are composed of many excitation levels of hydrogen 
molecules (H 2 ), and at high intensities the LW pho¬ 
tons can photo-dissociate H 2 . Because first stars are 
born with the help of H 2 -cooling, calculating LW back¬ 
ground is crucial in the study of first-star formation. 

The functional form of the RGF for L W background 
was first laid out by lAhn et all ( 200911 . The band- 
averaged and angle-averaged intensity Jlw(x, z) ob¬ 
served at comoving location x and redshift z due to 
a source at (x', Zg) in units of ergs“^ cm“^ Hz“® sr“^ 
is given by 


Note that once the box size, the grid resolution and 
the output times are determined, one can precalcu¬ 
late Gnii, j, k) for each set of observing time (r) and 
lookback time (t„-t„+i) and store all Gnip, q, f)’s re¬ 
gardless of the source distribution, which is a merit of 
the convolution theorem. The number of data files for 
Gnii, j, k^ is ^ A^^bserveA^lookback; where A^^^iserve 1® the 
number of observing redshifts and A^iookback is the num¬ 
ber of lookback times. Mookback can depend on the ob¬ 
serving redshift, depending on the change of tic over 
time. When the box size (L) is much smaller than the 
size of the effective horizon ivh), the number of unit 
sources repeated around the computation domain will 


Jlw(x, z) = 


1 

(47r)^ 


(Lt.(x', Zg)) (1 -k zf 

|x — x'l^ 1 + 2 s 


/mod (|x 


(14) 

where (Li,(x', z*)) = dihpiy)L,,ix', Zs)/(2.1eV) 

is the rest-frame, frequency(j^)-averaged source lumi¬ 
nosity of the emitted LW band photons at the source 
redshift Zs and hp is the Planck constant. 


/mod 



( 




116.29a 


0.68 


- 0.7 


if hzlll <97.39 

Q! — 

otherwise, 

(15) 
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Figure 2. Source distribution and LW background from different time slices, with the same configuration as in Figure [T] 
The top row shows the 2D maps of {Lv) (assumed to be fixed at and beyond the redshift specified below each panel until 
the adjacent future redshift) on a slice of thickness (114/256)//i Mpc, and the bottom row the “partial” Jlw’s on the same 
slice contributed from the look-back time ranges specified below each panel. The net Jlw observed at the observing redshift 
a = 27.90 is given by the summation of all the partial Jlw’s, which is shown in Figure . 


and 


a = 



0.27/ V 21 y 


(16) 


with the Hubble constant (in units of 100 km/s/Mpc) 
h and the present ratio of matter density to the critical 
density flm- Of course, the light-cone effect is implicit 
here: |x — x'| = t{z)—t{zs). The LW horizon is rpw = 
97.39a, beyond which no sources contribute to Jlw- 
Comparing Equation m to Equations © and dS]), we 
find that 


J-(|x 


c'l , T t') 


1 [l-fz(r)]^ 

(47r)^ |x — x'l^ 1 + ^s('r0 
x/mod(|x-x'|), (17) 


if we take the source term to be a sum of in¬ 
dividual source luminosities such that /(x', r') = 

Es Zs{t'))) Six' - x"). 

When visualized, the real-space RGF C/„(f, j, k) can 
give some insight on how Jpw is constructed. Figure [T] 
depicts RGFs for Jlw, which clearly show how the re¬ 
gion of influence is determined by the light-cone effect 
from a unit source at one corner of the boslfl. Follow¬ 
ing the steps described in Section [2?^ we have calcu¬ 
lated Jlw and used this to study the formation and 
suppres sion of the first stars d u ring the epoch of reion- 
ization ( Ahn et ^l2009ll2012ll . lAhn et ^ (|2912ll sim¬ 
ulated the process of cosmic reionization including the 
rise and fall of first stars, in a box of size 114//i Mpc, 
on which Figures [T]-|3] are based. Lifetime of first stars, 
which are believed to be as massive as M* > 100 Mq 


®The periodic boundary condition makes unit sources appear 
“around” every corner of a given box, as seen in Figure^ 


on av erage (e.g. see a recent result by iHirano et al.l 
12014 . is only about a few Myrs or less, and thus the 
life and death of these stars render the source func¬ 
tion to change rapidly in time and space. In Figure [U 
the rapid evolution of {L^ every 2 Myrs due to first 
stars is shown, together with the partial contributions 
Fn{p, q, r) to Jlw- ?’lw resides within the lookback 
time of 10 Myrs. One can see that contributions to 
Jlw becomes weaker and more homogeneous as look- 
back time increases. Nevertheless, contribution of those 
time slices with lookback time tiookback = [2 — 10] Myrs 
(2nd to 5th panels from left in Figures [T] and [2]) to Jlw 
is not completely negligible. 

How important is the light-cone effect in the LW 
background from first stars? To answer this question, 
we now calculate Jlw in an action-at-a-distance way, 
and then compare this to Jlw correctly calculated in¬ 
cluding the light-cone effect. For this, we now freeze the 
source distribution at the observing redshift, while at¬ 
taching identical boxes side by side for periodicity. Ba¬ 
sically, we now calculate this “instantaneously affected” 
Jlw as 


Jlw, instant (x, z') 


1 {L^{x', z)) {1 + zA 

(47r)^ jx —x'j^ 1 + 2s 

x/mod(|x-x'|), ( 18 ) 


where (LA is now frozen at the observing redshift z 
(e.g. the source term corresponding to the top-leftmost 
panel in Figure [5]), but we keep all other terms intact 
including the implicit condition |x — x'j = t{z) — t(zs), 
in order to quantify only the impact of the tempo¬ 
ral evolution of the source term. Figure [3] shows that 
the rapid evolution of source term indeed imprints its 
signature in Jlw- The fractional difference fdis = 
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Figure 3. Comparison of Jlw’s with the light-cone effect and without. Jlw under the correctly implemented light-cone effect 
(left), Jlw under an action at a distance (middle; see Equation 1 1811 and the fractional difference between these two fields 
(right) are plotted in the same 2D slice as in Figures [T] and [21 


(-^LW, instant “ ^Lw)/-^LW IS ubout -- 10%, which is quite 
significant. Therefore, astrophysical processes governed 
by radiation sources of short lifetime should be affected 
by the light-cone effect in general. 

4. Discussion 

We have formulated a generic method to numerically 
calculate radiation fields which are affected by the light- 
cone effect. This method requires discretizing both 
time and space coordinates, and is analogous to the 
usual mesh scheme used in N-body gravity simulations. 
The light-cone effect is realized by the retarded-time 
Greeen’s function, and its calculation becomes efficient 
when a periodic-box condition is used. Once the RGFs 
are calculated, the process of obtaining a radiation field 
of interest is performed by FFT, which greatly expe¬ 
dites computation compared to direct summation over 
all sources inside and outside the simulation box. 

The relevance of this method is clear in astrophysical 
problems occurring on cosmological scales. Once the 
size of the simulation box is chosen and the functional 
form of the RGF is known, one can follow the algorithm 
in Section as in the example of LW background in 
Section |31 Similar large-scale astrophysical processes 
which require an accurate account of the light-cone ef¬ 
fect are the reionization of baryons by UV photons and 
the reionization and heating of gas by X-ray photons. 
It becomes more important at very high redshift when 
most astrophysical sources are first stars and their by¬ 
products, because the lifetime of first stars are believed 
to be very short and thus the source function evolves 
rapidly in time (Section [3]). 

There are, of course, limitations to this method. The 
method is based on the fact that the RGF is a function 
only of distance to the source but not of direction to the 
source. In reality, radiative transfer is usually affected 
by the actual path of the ray in the form of the locally- 
varying optical depth. For example, formation of an H 
II bubble from a source does not occur in a spherically 
symmetric way, even though light travels at the speed 
of light. This is due to the inhomogeneous distribu¬ 


tion of hydrogen density (and also the optical depth), 
which then regulates the H-ionizing background to be 
dependent on the direction of the ray. For H-ionizing 
background, our method can then be used only after the 
IGM is fully ionized and the optical depth becomes neg¬ 
ligible everywhere. In case of LW background, we could 
safely use the method because the modulation factor 
/mod in Equation (iTSl) . derived from the cosmic red- 
shifting of LW photons and the frequency of hydrogen 
Lyman resonance lines, is valid even when the IGM is 
ionized, because simulations of cosmic reionization find 
that the IGM leaves trace amount of neutral hydrog en 
atoms due to recombination (e.g. Illiev et all (l200^ ). 
In this case, very-high opacity of Lyman resonance lines 
guarantees that once LW band photons are redshifted 
into one of these lines, even a trace amount of hydro¬ 
gen atoms can successfully scatter off those photons, 
and then /mod becomes the only modulation factor in 
addition to the usual geometrical dilution factor (Equa¬ 
tions |TT] and [T7]). Another limitation is that the method 
cannot treat the gravitational lensing effect, which de¬ 
pends on the actual mass distribution around each ray’s 
path. If a problem of interest is found to be strongly 
affected by the lensing effect, one should resort back to 
the ray-tracing method. In addition, if source locations 
are resolved beyond the mesh resolution, one can adap¬ 
tively choose to use the particle-mesh (PM) scheme to 
calculate the near-source influence by direct summa¬ 
tion over nearby-source contributions while calculating 
far-source influence by our method. 
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